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We develop an empirical potential for silicon which represents a considerable improvement over 
existing models in describing local bonding for bulk defects and disordered phases. The model 
consists of two- and three-body interactions with theoretically motivated functional forms that 
capture chemical and physical trends as explained in a companion paper. The numerical parameters 
in the functional form are obtained by fitting to a set of ab initio results from quantum mechanical 
calculations based on density functional theory in the local density approximation, which include 
various bulk phases and defect structures. We test the potential by applying it to the relaxation 
of point defects, core properties of partial dislocations and the structure of disordered phases, none 
of which are included in the fitting procedure. For dislocations, our model makes predictions in 
excellent agreement with ab initio and tight-binding calculations. It is the only potential known 
to describe both the 30°- and 90°-partial dislocations in the glide set {111}. The structural and 
thermodynamic properties of the liquid and amorphous phases are also in good agreement with 
experimental and ab initio results. Our potential is the first capable of simulating a quench directly 
from the liquid to the amorphous phase, and the resulting amorphous structure is more realistic than 
with existing empirical preparation methods. These advances in transferability come with no extra 
computational cost, since force evaluation with our model is faster than with the popular potential 
of Stillinger- Weber, thus allowing reliable atomistic simulations of very large atomic systems. 



I. INTRODUCTION 



Silicon is one of the most intensely investigated mate- 
rials. Interest in silicon is mostly motivated by its great 
technological importance, but it is also regarded as the 
prototypical covalent solid, on which theoretical concepts 
about covalent bonding can be tested, and new ideas can 
be explored. The nature of covalent bonding often makes 
the description of complicated phenomena difficult: the 
covalent bond can only be described properly by quan- 
tum theory, while many interesting phenomena involve 
large numbers of atoms (in the range 10^-lOE^, which 
quantum mechanical approaches cannot handlea. In this 
sense, silicon represents an ideal candidate for modeling 
by an empirical interatomic potential: the development 
of a potential will at once put to a rigorous test our un- 
derstanding of the physics of covalent bonding and, if 
successful, will enable the simulation of important com- 
plex processes that involve large number of atoms. 

Many attempts to construct an empirical model for i: 
teratomic interactions in Si have already been reportC' 
Of all these, models, the Stillinger and Weber (SW|3 and 
the TersoficTO potentials are the most widely useda. The 
SW potential consists of two- and three-body terms and 
was fitted to experimental propertieSj-pf the diamond cu- 
bic (DC) and molten phase of silicon0_.1t has been u! 
for exarpBle, to study lattice dynamicsQ, point de 
surfacegl3, and the liquid and amorphous phasesP'''^^' 



The Terseff potential (its Ahree versions usually referred 
to as TlQ, T20, and T3H) consists of many-body in- 
teractions included in a bond order term and was fit- 
ted to ab initio results for several Si polytypes. It has 
been usedrto study lattice |disnamicsQ, thermomechanical 
propertiest3, Jioipt, defectaJ^tl, and the liquid and amor- 
phous phaseaj'tja. 

Although the SW and Tersoff functional forms have 
enough flexibility to describe a number of different config- 
urations, their transferability to a wide class of structures 
remains in question. Several models have attempted to 
improve the description of configurations far from equi- 
librium and far from the database used to construct the 
potential, by changing the functional form, using higher 
order (up to five-body) expansion terms, increasing the 
number of fitting parameters .-(up to 36), or expanding 
the set of fitted structure3l3~E3. Despite the increased 
complexity, these models have been able to improve the 
description of local configurations only selectively, and 
often at a considerable increase in computational cost. 
This suggests that a simple extension to more elaborate 
functional forms or larger databases does not necessarily 
provide better description of covalent bonding. Consid- 
ering the lack of transferability of existing models, it is of 
interest to develop a model for silicon with the following 
ingredients: superior description of local structures and 
computationally efficient evaluations of the energy and 
interatomic forces. 
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In this article we present a new empirical potential for 
silicon using a theoretically motivated functional, form 
which emphasizes chemical and physical trendsEll, and 
which is determined by fitting to a fairly small ab initio 
database with only 13 parameters. This potential rep- 
resents a considerable improvement over existing models 
in describing local structures and extended defects. It 
provides a good description of point defects in the bulk, 
the concerted exchange path for self-diffusion, and elastic 
properties of bulk silicon. It also predicts core structures 
of partial dislocations in the glide set {111} in excellent 
agreement with ab initio results. Disordered structures 
and phase transitions are also well-described, particularly 
the amorphous phase, which is better modeled by dynam- 
ical simulations using our potential than by any empirical 
preparation method. 

The article is organized as follows: The functional form 
is presented in detail in section ||. The fitting of the 
model is discussed in section |l|, along with tests of trans- 
ferability for bulk crystal structures, defects and activa- 
tion complexes. The fitted potential is then used to cal- 
culate core properties of physically relevant dislocations 
in section |l^ and the structure of disordered phases in 
section ^ 



II. FUNCTIONAL FORM 

Here we describe the functional form of the new po- 
tential, which we call the Environment Dependent Inter- 
atomic Potential (EDIP) for bulk silicon, and refer the 
reader to the companion papeicj for theoretical justifi- 
cation of all the terms. The energy of a configuration 
{Ri} is expressed as a sum over single- atom energies, 
E = Ei, each containing two- and three-body terms, 
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where V2{Rij , Zi) is an interaction between atoms i and 
j representing pairwise bonds, and V^lRij, Rik, Zi) is the 
interaction between atoms i, j ^-nd k centered at atom 
i representing angular forces. Both types of interactions 
depend on the local environment of atom i through its 
effective coordination number, defined by. 
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where f(Rim) is a cutoff function that measures the con- 
tribution of neighbor m to the coordination of atom i 
in terms of the separation Rim- The neighbor function, 
depicted in Fig. |l|, is unity for r < c, with a gentle drop 
from 1 to between c and a, and is for r > a: 




if r < c 
(t^) ifc<r<a 
if r > a 



(3) 



■ A neighbor of atom i at a distance 
r < c is considered a full neighbor, while the neighbors 
between c and a give only a partial contribution to Zi. 
The cutoffs are constrained to reproduce the coordina- 
tions of important crystal structures, e. g. Zi = A for 
the diamond lattice. 

The two-body term includes repulsive and attractive 
interactions, 



V2{r, Z)^A 



exp 



(4) 



which go to zero at the cutoff r = a with all deriva- 
tives continuous. Although this choice is quite similar 
to the SW two-body term (and indeed reduces to it ex- 
actly for small distortions of the diamond lattice), the 
bond strength adapts to changes in the local atomic 
environment. The coordination dependence introduces 
an asymmetry, V2{Rij, Zi) ^ V2{Rji, Zj), similar to the 
Tersoff potential in that the strength of the attractive 
force is controlled by a bond order function p{Z) that 
depends on the local coordination. This dependence is 
motivated by theoretical calculations which have demon- 
strated the weakening of the attractive interaction and 
the corresponding increase in bond length for increasing 
coordinationEll a. These studies reveal a shoulder in the 
function p{Z) at the ideal coordination Zq — 4, where the 
transition from covalent {Z < Zq) to metallic {Z > Zq) 
bonding occurs. This theoretical dependence can be ac- 
curately captured by a Gaussian functiorciJ, 



piZ) = c 



(5) 



Figure || shows the V2 (r, Z) term for several coordina- 
tions, compared to the two-body part of the SW potential 
V2^^{r). Note that the fitted V2{r, Z) resembles closely 
the inverted ab initio pair potentials for silicon crystals 
with the sam£_.coordinationg23, a feature built into the 
function formed. 

The three-body term contains radial and angular fac- 
tors. 



Vz{Rij,Rik,Zi) = g{Rij)g{Rik)h{lijk, Zi), 
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where lijk = cos 



Rij • Rik/RijRik- The radial 



function is chosen to have the SW form, 



g{r) = exp 
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and goes to zero smoothly at the cutoff distance a. The 
angular function h{l, Z) is strongly dependent on the lo- 
cal coordination through two functions t{Z) and w{Z) 
that control the equilibrium angle and the strength of 
the interaction, respectively. TheoretiGaJ considerations 
lead us to postulate the following formEil: 
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where H{x) is a generic function satisfying the con- 
straints, H{x) > 0, H{0) = 0, iJ'(O) = and H'' (0) > 0. 
Our specific choice is 

h{i, z) = \ [(i ~ e-Q(m+r[z)?^ ^ vQim + T{z)f 

(9) 

We make the choice w{Z)^'^ = Q{Z) = Qae^f^^ to con- 
trol the strength of angular forces as function of coordi- 
nation. The three-body angular function becomes flat- 
ter (and hence angular forces become weaker) as coor- 
dination increases, representing the transition from co- 
valent bonding to metallic bonding. The angular func- 
tion /i(Z, Z) has two major contributions. The first, 
Hi{x) oc 1 — e~^ , is symmetric about the minimum, 
becoming flat for small angles. This choice of shape for 
the angular function was .also used by Mistriotis, Flyntza- 
nis and Farantos (MFF)Ea, but due to its environment- 
dependence our angular function is fundamentally dif- 
ferent. In a preliminary fitted version-of EDIP, wc in- 
cluded only this flat, symmetric termed, but we found 
that it is not suitable for several structures (particu- 
larly the amorphous phase). Indeed, a more asymmet- 
rical angular function is suggested by approsiiaations of 
quantum-mechanical (tight-binding) modelsc3"E3 aud ex- 
act inversions of ab initio cohesive energy curve^J, so 
the present form of EDIP also contains a second term, 
H2{x) oc x"^ (identical in shape to the SW form but con- 
taining environment-dependence) , which gives a stronger 
interaction for small angles. 

The function t{Z) = ~1(){Z) ^ - cos(6'o(Z)) controls 
the equilibrium angle 9q{Z) of the three-body interaction 
as a function of coordination. This feature of the poten- 
tial makes it possible to model the. proper hybridization of 
atoms in different environment£j'Lj: If a silicon atom is 
three- or four-fold coordinated, it will prefer to form sp^ 
or sp^ hybrid bonds with equilibrium angles 0o(3) = 120° 
and 6*0(4) = 109.471°, respectively. For coordinations 2 
and 6, we take r(2) = t(6) = or 6*0(2) = 6lo(6) = 90°. 
For two-fold coordination, this choice describes the pref- 
erence for bonding along two orthogonal p-states with 
the low-energy nonbonding s state fully occupied. For 
six- fold coordination, the choice 9q{Q) = 90° reflects 
the p character of the bonds (which are also metallic, 
as discussed below). We construct t{Z) to interpolate 
smoothly between the special points {Z = 2, 3, 4, 6) with 
the following form. 



t{Z) = ui -I- ^2(^36' 



(10) 



with the parameters chosen as u\ — —0.165799 , ui = 
32.557, u-i = 0.286198, and -U4 = 0.66, resulting in 
the curve shown in Fig. |^, which is consistent witli 
the results of quantum-mechanical approximations^. 
Note that these parameters are theoretically determined 



and are not allowed to vary in the fitting described 
in the next section. Figure |4| shows the three-body 
term Vs{Rij , Rik, Zi) for three atoms at a distance 2.35 
A, and for different coordinations. We also compare 
the three-body term to the SW three-body potential 
Vi^{Rij,R,k)- The SW angular form hsw{Ojik) = 
^Swicos{9jik) + 1/3)^ penalizes the configurations with 
angles smaller than 90° with a large positive contribution. 
In contrast, the angular function of our model potential 
gives a considerably weaker interaction at small angles. 

In summary, this implementation of EDIP for bulk sil- 
icon has 13 adjustable parameters: A, B, p, (3, a, a, c, 
A, 77, 7, Qo, A* ^-iid a. It also has continuous first and 
second derivatives with respect to the atomic position 
vectors. The functional form already contains consider- 
able information about chemical bonding in bulk silicon 
taken directly from theoretical studies, mostly of ideal 
crystal structures. The (relatively few) adjustable pa- 
rameters provide the necessary freedom to extrapolate 
these bonding dependences for defect structures strictly 
outside the theoretical input, as described below. 

We close this section by addressing the crucial issue 
of computational efficiency (which motivates the use of 
empirical methods in the first place). The environment 
terms in the two- and three-body iterations require extra 
loops in the force calculation. In the case of the three- 
body interaction, it requires a four-body loop, that would 
make a force evaluation more expensive than with the 
SW potential. However, the four-body loop needs to be 
performed only for those neighbors I of atom i in which 
df{Rii)/dRii 0. This happens only when atoms are 
in the range c < r < a, i.e. only for a small number 
of the neighbors. Therefore, one force evaluation using 
our model is approximately as efficient as one using the 
SW potential, showing that increased sophistication and 
realism can be achieved with insignificant computational 
overhead. In fact, since the fitted cutoff distance (given 
in the next section) is smaller than the corresponding 
SW cutoff, computing forces with our model is typically 
faster than with the SW potential. For example, in liquid 
phase simulations on a Silicon Graphics R-10000 proces- 
sor it takes 30 ^sec/atom to evaluate forces using EDIP, 
compared with 50 //sec/atom using the SW potential. 



III. FITTING AND TESTS 

In order to determine the values of the adjustable pa- 
rameters, wc fit to a database that includes ab initio 
resulta23, based on density functional theory in the local 
density approximation (DFT/LDA), for bulk properties 
(cohesive energy and lattice constant of the DC struc- 
ture), selected values along the unrelaxed concerted ex- 
change (CE) pathEll for self-diffusion, some formation en- 
ergies of unrelaxed point defects (vacancy and interstitial 
at the tetrahedral and hexagonal configurations )E2rEJ, a 
few key values in the generalized stacking fault (GSF) en- 
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ergy surfaced, and the experimental elastic constantscS. 
This rich set of configurations contains many of the im- 
portant local structures found in condensed phases and 
bulk defects and thus improves the transferability of the 
model. The database does not include other high sym- 
metry configurations, such as SC, BCC or FCC (as many 
of the existing empirical potentials have done) , although 
the DC structure is required to be the lowest in energy. 
These hypothetical metallic structures have large enough 
energies compared to covalent ones to be considered ir- 
relevant. Incidentally, fitting the model simultaneously 
to such a wide set of configurations and properties repre- 
sents a considerable computational challenge. The fitting 
is accomplished using a least-squares approach, with each 
configuration in the database weighted appropriately. All 
parameter values are allowed to vary at once through 
a simulated annealing process. Theoretical estimates of 
the parameters were also used to restrict the. range of 
parameter-space that needed to be exploredEJ. Table | 
provides the best set of parameters we obtained. 

The calculated energies and properties of several bulk 
structures as obtained from our potential and from ab ini- 
tio and other empirical potential calculations are given 
in Table ||. This compilation of results includes the 
DC structure which is the ground state of Si, and sev- 
eral other high-coordination bulk structures. Although 
the latter structures were not included in the fitting 
database, our model provides a reasonable description 
of their energies and a good description of equilibrium 
lattice constants. Structures such as /3-tin and BCT5lD 
are also of interest because they have low-energy and rel- 
atively low-coordination (5 for BCT5 and 6 for /3-tin)El. 
Experimentally, the DC phase transforms into the /9-tin 
structure under pressure. For the /3-tin structure, our 
model predicts a cohesive energy per atom higher than 
that of DC by AE — 0.67 eV, and a lattice parameter 
Go = 4.76 A, as compared to ab initio results of 0.21 
eV and 4.73 A respectively. For the BCT5 structure, 
our model predicts AE = 0.26 eV and Oo = 3.36 A, 
as compared to ab initio results of 0.23 eV and 3.32 A, 
respectively. 

Most existing empirical potentials give a poor (or 
marginally acceptable) description of elastic properties 
of the DC crystal, which directly affects the description 
of the crystal deformation. In the fitting database we in- 
cluded the three independent elasti c co nstants Cn, C12 
and C44 from experiment. Table III compares elastic 
constants,, obtained with the homogeneous deformation 
methodE3, as given by our model with the results from 
other empirical potentials and with experimental results. 
The shear constant C44 is particularly, important for the 
description of long range interactionso. Although C44 is 
underestimated by most empirical models, ours is rea- 
sonably close to the experimental value. Our potential 
also almost perfectly reproduces the ab initio vahie of 
the shear modulus C4^without internal relaxationLJ, as 
anticipated by theorjO. Table III also includes other 



elastic properties, such as the second shear constant. 



Cii — C12, and the Cauchy discrepancy, AC = Cig-j C44, 
both important for determining crystal stabilityc3. By 
predicting a negative Cauchy discrepancy, our potential 
offers a qualitative improvement over most other exist^ 
ing potentialsci and several tight-binding (TB) modelsEll 
which give positive values. It also provides a quantitative 
improvement, over other TB models that get the correct 
sign of ACll3. In summary, our potential gives elastic 
constants in excellent agreement with experimental and 
ab initio results. As we have shown in Ref.Ej, accurate 
elastic behavior is not simply the result of good fitting, 
but is rather a built-in feature of our functional form. 

Point defects in the DC crystal involve large atomic re- 
laxations and rebonding, thus representing the first test 
for the transferability of our model in describing local 
structures away from equilibrium. Our fitting database 
included the unrelaxed structures of the vacancy (V) and 
the interstitialJn tctrahcdral (It) and hexagonal (Ih) 
configurationsEB. Table IV shows the formation energy 
for the unrelaxed and relaxed structures of V, It aiuil#. 
as obtained using oiit,model, compared to ak-.initi(x3'l3 
and TB calculationalll, and from calculations^ using the 
SW and Tersoff (T2 and T3) potentials. The relaxed 
structures are computed lising an energy minimization 
conjugate-gradient methods. Although the SW and Ter- 
soff potentials give a reasonable description of relaxed 
structures, they clearly fail in describing the energy re- 
leased upon relaxation. In our model, on the other hand, 
the relaxation energies are much lower and in reasonably 
good agreement with ab initio calculations. We empha- 
size that none of the relaxation energies or structures 
were used in the fitting database. The < 110 > split 
interstitial (the lowest energy interstitial configuration) 
is well described by our model, in spite of also not be- 
ing included in the fitting. The formation energy for the 
relaxed structure of the < 110 > split interstitial is 3.35 
eV in our model, compared to 3.30 eV from ab initio 
calculationala, while-.the SW potential gives 4.68 eV. 

The CE processEil has been ideatified as a possible 
mechanism for self-diffusion in Sil-3. Most of the em- 
pirical potentials provide only a fair description of this 
important and complicated sequence of configuratiDnalJ. 
Fig. 1^ shows the energy for the unrelaxed CE pathO from 
calculations based on DFT/LDAE2I, the present model 
and the SW and Tersoff potentials. The results using our 
model agree reasonably well with those from DFT/LDA 
calculations, and are considerably better than those us- 
ing the SW or Tersoff potentials. The activation energy 
for CE (an important quantity that enters in the calcu- 
lation of diffusion rates) obtained from our model is 5.41 
eV, compared to 5.47 eV from ab initio calculations, 7.90 
eV from the SW potential, and 6.50 eV from the Tersoff 
potential. The energy of the relaxed structure, which 
was not included in the fitting database, is 4.82 eV using 
our model, in good agreement with the ab initio result of 
4.60 eV. 

The fitting database also included selected unrelaxed 
configurations of the generalized stacking fault (GSF) en- 
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ergies. We considered threCp^oints for the the ghde set 
and three for the shuffle selE3 of the {111} ghde plane. 
Table ^ shows the unstable stacking fault energy for the 
unrelaxed glide and shuffle {IfU} planes as obtained from 
calculations using DFT/LDAE3, SW, and the present in- 
teratomic potential. Our model gives good agreement 
with DFT/LDA calculations for the < 112 > ghde set, 
but underestimates the energy for the < 110 > glide set. 
Since our potential is short-ranged, it gives zero energy 
for the stable stacking fault, while the experimental value 
is 0.006 eV/A2. The SW and Tersoff models, which are 
also short-ranged, give the same zjeKO energy for the sta- 
ble stacking fault. The TB modeled, on the other hand, 
gives 0.005 eV/A^, in agreement with DFT/LDA results. 
Nevertheless, given that the accuracy of empirical models 
is rarely better than a few tenths of an eV per atom, our 
vanishing stable stacking fault energy is not particularly 
problematic. 

In summary, our potential provides an excellent de- 
scription of configurations near equilibrium as well as a 
wide range of point defects in the DC structure. Al- 
though a number of these properties were explicitly fit, 
the superior description has been achieved with a small 
number of adjustable parameters that is greatly exceeded 
by the number of degrees of freedom inherent in these 
configurations. Thus, we suggest that it is our physi- 
cally appropriate functional form rather than a flexible 
fitting strategy that has led to the improved description, 
a conclusion that is further supported by our results for 
extended defects and disordered phases discussed next. 



IV. DISLOCATIONS 

A stringent test of the transferability of our model to 
local structures is obtained by calculating core properties 
of dislocations, about which no informatiQii|jii|as included 
in the fitting database_-Several ab mitidfJcj and tight- 
binding calculationsE£lL2l have been performed for dislo- 
cations in silicon. Although such calculations are feasible 
only for small systems containing of order few hundred 
atoms, they provide important information,|-such as the 
core structure of the 90°-partial dislopationcZl, the kink 
structure iuptlie 30°-partial dislocationo, and dislocation 
interactionsEJ. However, full simulations of dislocation 
dynamics and its effect on the mechanical properties of 
the solid require much larger cells owing to the long range 
interaction of the stress fields, and therefore can only be 
performed using methods that are computationally less 
expensive. Empirical models have beep-Jised to study 
several aspects of dislocations in silicoroEj'LJ. No such 
model has prpven reliable enough to provide a description 
of long rangeEj as well as core propertiesLj of dislocations 
at the same time. The flaw in describing long range ijiter.- 
actions is due to a poor description of elastic forcesEJo, 
while the flaw in describing core ppoperties is due to a 
poor description of local structured^. Existing models 



do not give the correct structure for reconstructed cores: 
for example, the SW potential gives reconstruction only 
for the 30°-partial dislocation, and the Tersoff potential 
gives reconstruction only for the 90°-partial. 

In the present study, all dislocation struc±ures are com- 
puted-using energy minimization methodscj at constant 
stressE3 for a system of 3600 atoms, and the cell bound- 
aries lie in the [112], [111], and [TlO] directions. Fig. ^ 
shows (a) the unreconstructed and (b) the reconstructed 
core structure of a 90°-partial dislocation. The recon- 
struction energy is deflned as the energy gain (per unit 
length) for the system to go from unreconstructed to re- 
constructed configuration. Table ^ compares the recon- 
struction energy in units of eV/6 (where b is the unit 
length of the dislocation, b — 3.84A), as given by our 
model and by SW, Tersoff, TB, and ab initio calcula- 
tions. Configuration (b) is neither stable nor metastable 
for the SW potential, i.e. the SW model does not sup- 
port reconstruction for the 90°-partial dislocations. The 
present model predicts configuration (b) as the lowest 
in energy, while configuration (a) is metastable. This 
model gives reconstruction energy of 0.84 eV/&, in excel- 
lent agreement with the ab initio value of 0.87 eV/ffiZl. 
In our calculation, the recpestructed rbonds are stretched 
by 2.1 %, while ab initial! and TBE2I calculations give 
bonds stretched by 2.6 % and 3.0 % respectively. 

Fig. ^ shows (a) the unreconstructed and (b) the re- 
constructed core structures of a 30°-partial dislocation. 
The results -for the reconstruction energy from ab initio 
calculationsEj, and from empirical potential calculations 
using the SW, Tersoff and our model, are presented in ta- 
ble Kg The Tersoff potential gives negative reconstruc- 
tion energy, so that the unreconstructed configuration 
(a) would be the more stableEj, contrary to experimen- 
tal and ab initio results. Although the SW model gives 
the correct reconstructed configuration, the reconstruc- 
tion energy is twice as large as the ab initio result. Our 
model, on the other hand, gives the reconstruction energy 
in very good agreement with the ab initio calculation. We 
find that the reconstructed bonds in the dislocatiau. core 
are stretched by 3.6%. The ab initio calculationsEj give 
bond stretching of 4.2% at most. 

It is important to mention that we found a few 
metastable partially reconstructed configurations be- 
tween the unreconstructed and reconstructed configura- 
tions for both the 90°- and 30°-partial dislocations. Such 
configurations are artifacts of our description of the local 
environment during changes ip,coordination, and proba- 
bly have no physical meaningH. 

A defect in the core of a recor^ructed dislocation is 
cafled an antiphase defect (APD)e3. Fig. | shows APD 
configurations in (a) a 30°-partial dislocation and (b) a 
90°-partial dislocation. Table ^ gives the correspond- 
ing APD formation energies. Since the SW model does 
not produce a reconstructed configuration for 90°-partial 
dislocations, there is no stable APD configuration for 
this case. For a 30°-partial dislocation, the SW model 
gives an APD formation energy much larger than ab ini- 
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tio calculationsE3. The Tersoff potential gives a negative 
value for the APD formation energy of the 30°-partial 
and a considerably smaller energy for the APD in a 90°- 
partial. Our model gives the APD formation energy in 
good agreement with ab initio calculations for the 30°- 
partial dislocation. For the 90°-partial dislocation, to 
our knowledge there is no ab initio APD energy calcu- 
lation available, but our APpQ energy is somewhat low 
compared to TB calculation£3. 

From the above comparisons we have established that 
our empirical model is the first to give a full description 
of core properties of both the 90° and 30°-partial dislo- 
cations in silicon with reasonable accuracy. Our model 
predicts reasonable reconstruction energies and provides 
a good description of local deformations in the disloca- 
tion cores. These results demonstrate remarkable trans- 
ferability since the local atomic configurations present in 
dislocation cores are quite different from the structures 
included in the fitting database. 



V. DISORDERED STRUCTURES 

Another stringent test of the transferability of our 
model for bulk material is the calculation of structural 
properties of the liquid and amorphous phases. Such 
disordered structures contain a rich set of local bond- 
ing states from covalent (amorphous) to metallic (liq- 
uid) about which no information was included in our fit- 
ting procedure. Existing environment-independent po- 
tentials have had considerable difficulty in simultane- 
ously desc.rihicLg the crystalline, liquid and amorphous 
phasea3'OE2L3. In the preceding sections we have 
demonstrated an improved description of the crystalline 
solid and its defects, so we now turn our attention to 
whether our environment-dependence can extrapolate 
these successes to the liquid and amorphous phases. 

Liquid phase. The SW potential has been shown 
to rejpmduce the pair correlation function g(r) of the 
liquidBEJ. It also predicts the melting temperature Tm to 
within a few hundred degrees of the experimental value 
of i685 kI (although it was explicitly fit to reproduce 
T„p). In spite of these successes, the SW poi£iitial has 
difficulty in describing the structure of liquidEj, i.e. it 
does not reproduce the ab initio bond-pajigle distribution, 
overly favoring angles near tetrahedraES. The Tersoff po- 
tentials, on the other hand, predicl||a g(r) that favors the 
unphysical four- fold coordinationOu, and the fwly ver- 
sion predicting reasonable liquid structure is T3Qo. The 
melting temperature, howesier, is greatly overestimated 
by T3 at around 3000 KBH. The first-neighbor bond 
angles and coordination statistics are not well described 
by T3, although the statistics are improved by using an 
(arbitrarily) longer, coordination cutoff beyond the first 
minimum of g(r)llj. We emphasize that the T2 poten- 
tial, the most-, successful parameterization of the Tersoff 
model overalB, cannot decribe the liquid phased. 



We prepared a 1728-atom liquid sample with the 
present potential dXT — 1800 K and zero pressure using 
standard simulation techniques, although we used a con- 
siderably lon ger t ime and larger system size than in pre- 
vious studiesEilCa. The structural properties of the model 

results 



(which 



liquid are shown in Fig. ^ and compared with tl 
of a 64-atom ab initio molecular dynamics studyl^ 
are similiar to recent results with 343 atoms including 
electron spin effects and gradient correctionsEj) . The pair 
correlation function g(r) shows excessive short-range or- 
der with our potential, as evidenced by the overly sharp 
first neighbor peak containing around 4.5 first neighbors, 
smaller than the experimental value of 6.4. This is con- 
sistent with the fact that in our model the density of 
the liquid is somewhat smaller than that of the solid, 
while in—reality silicon becomes 10% more dense upon 
meltingEfl. Although these features are unphysical, the 
present model offers a qualitative improvement in the 
bond angle distribution function 33(0, r^), which gives 
the (normalized) distribution of angles between pairs of 
bonds shorter than r^, the first minimum of g{r). As 
shown in Fig. ^, our potential predicts the auxiliary 
maximum a,i 9 = 60°, although the primary maximum 
is shifted toward the tetrahedral angle away from the 
ab initio most probable angle oi = 90°. The present 
model is the first to capture the biraOjdal shape of the 
first-neighbor bond angle distributionE^. 

The thermodynamic properties of the melting transi- 
tion (aside from the change in density) are reasonably 
well-described by our model, in spite of its not having 
been fit to reproduce any such quantities. The bulk melt- 
ing point Tjn predicted by our model is 1370 ± 20 K, 
which is 20% below the experimental value. The melting 
point was measured for a finite sample with (100)2 x 1 
surfaces, heated from 300 K to 1500 K in 2 ns (over 10 
million time steps). A bulk solid with periodic bound- 
ary conditions superheats and melts around 2200 K. The 
latent heat of melting is 37.8 kJ/mol, in reasonable agre- 
ment with the experimental valuCrOf 50.7 kJ/mol, closer 
than the SW value of 31.4 kJ/molB. 

In summary, although the liquid has some unphysical 
features with our potential, it offers some improvements, 
particularly in describing bond angles. It is important 
to emphasize that reasonable liquid properties are pre- 
dicted by our model without any explicit fitting to the 
liquid phase; in contrast, the only two potentials reported 
to give an adequate description of the liquid, SW jaad 
MFF, were each fit to reproduce the melting pointBEa. 
With our model, the reduced density and excess of co- 
valent bonds may be artifacts of the short cutoff of our 
potential, which is appropriate for the covalently-bonded 
structures used in the fitting, but is perhaps too short 
to reppxiduce overcoordination in metallic phases like the 
liquidEl. 

Amorphous phase. Experimentally, amorphous silicon 
is known to form a random tetrahedral network, with 
long-range dispjder and short-range order similar to that 
of the crystalO'Ea. Ab initio molecular dynamics sim- 
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ulations of quenching a 64pa.tom liquid predict almost 
97% four-fold coordinatior£3. An empirical potential 
would be invaluable in exploring larger system sizes and 
longer relaxation times than are feasible from first prin- 
ciples, but unfortunately no existing potential is capable 
of quenching directly from the liquid to the amorphous 
phase. Instead, empirical model liquids typically trans- 
form into glassy phases |Uppa|-coaling, characterized by 
frozen-in liquid structuralJilj'E2rE2l. Alternatively, in the 
case of the T2 potentiaL-jquenching results in a reasonable 
amorphous structuraaO, but the original liquid phase is 
not riealistic and already contains excessive tetraherdal 
orderO. Therefore, it has been impossible to simulate 
an experimentally relevantjiath to the amorphous struc- 
ture [e.g. laser quenchingjlJ) , and artificial preparation 
methods have been required to create large 
phous structures with empirical potentialstirll^ 
For example, the so-called indirect SW amorphous (ISW) 
structure is created by increasing the strength— of the 
SW three-body interaction during the quenchlljil3. The 
ISW structure is stable with the unaltered SW potential, 
but since it has 81% four-fold coordinated atomaiil, it 
only beais, a weak resemblance to the real amorphous 



structuraHj. Abandoning molecular dynamics, an im- 
proved amorphous phase with close to 87% four-fold co- 
ordination can be generated using the band-switching 
algorithm of Wooten, Winer and WeairolJ, but such a 
method does not permit accurate simulation of atomic 
motion. A more realistic, large-scale amorphous struc- 
ture can be prepared using a hybrid of these methods: 
Holender and Morgan create an amorphous structure of 
over 10^ atoms with almost 94% four-fold coordination 
by patching together a number of smaller WWW struc- 
tures, thermally treating the sample at high tempera- 
tures and then relaxing it using molecular dynamics with 
the SW potential modified for stronger angular forcesE3. 
As these authors emphasize, however, this preparation 
method (which we call the HM2 model) was designed 
by trial-and-error to fit the experimental structure fac- 
tor and bears no resemblance to the real experimental 
generation of a-Si. They also report that if the SW po- 
tential is not modified, their preparation method results 
in a structure (which we call the HMl model) with only 
74% four-fold coordination. 

Remarkably, the present model predicts a quench di- 
rectly from the liquid into a high-quality amorphous 
structure. The phase transition is quite robust, since it 
occurs even with fast cooling rates. For example, quench- 
ing at -300 K/ps leads to a reasonable structure with 84% 
four-fold coordination. At much slower quench rates of 
-1 K/ps, an improved structure of 1728 atoms at T = 300 
K and zero pressure is produced with almost 95% four- 
fold coordination. The excess enthalpy of the amorphous 
phase compared with the crystal is 0.22 ay./atom, closer 
to experimental values < 0.19 eV/atonJl3 than the ah 
initio value 0.28 eV/atom (probably due to the con- 
strained voktfne and small system size used in the ab 
initio studyM). 



The coordination statistics of the am orph ous phase ob- 
tained with ouK-opodel, given in Table VII, are closer to 
ah initio resultsE^I than with most of the empirical models 
described above. The HM2 model provides a compara- 
ble description, but we stress that its preparation proce- 
dure is unphysical and that the modification of the SW 
potential necessary to achieve the improved descri ption 
(see the difference between HMl and HM2 in Table |vi| ) 
degrades many important properties, such as elastic con- 
stants, defect formation energies and the melting point. 
Since realistic preparation methods and dynamics have 
not been achieved with interatomic potentials, in the fol- 
lowing we compare our results only with experiments and 
ab initio simulations. 



The pair correlation function shown in Fife. 10 (a) 
is in good agreement with ab initio resultsE3. More- 
over, Fig. ^ shows that the radial distribution function 
t{r) — A'Kprg{r) is in excellent agreement with the ii^ults 
of neutron scattering experiments by Kugler et. aliEB (us- 
ing their experimental density p = 0.054 atoms/ A'^ for 
comparison). The persistence of intermediate-range or- 
der up to 10 Acaptured by our model as in experiment 
is a strength of the empirical approach, since this dis- 
tance is roughly the size ofr-the periodic simulation box 
used in the ab initio studiei^. Given the limited resolu- 
tion of the experimental data, especially at small r (large 
q in the structure factor), the sharper first three peaks 
with our model may be interpreted as refinements of the 
experimental results. In Table Vlll we summarize a de- 



tailed comparison of features of a-Si as obtained with our 
model and from ab-initio results, against those revealed 
by experiment. Overall the agreement between experi- 
ment is very satisfactory, with the results of the present 
model somewhat closer to experimental values than ab 
initio results as in the case of the enthalpy (AHa-c) and 
the bond-length (cr,i) and bond-angle (ag) deviations. 

The bond angle distribution 33(6*, r^) shown in Fig. 
(b) is narrowly peaked just below the tetrahedral an- 
gle, and also reproduces the small, weU^eparated peak 
at 60° observed in ab initio simulationsEj (unlike in pre- 
vious empirical models). The average angle is 108.7° and 
standard deviatioLLL,13.6° in close agreement with the ex- 
perimental valuesB of 108.6 ±0.2° and 10.2 ±0.8° and ab 
initio valued^ of 108.3° and 15.5°, respectively. Notice 
that the peaks in both g{r) and g^{9^rm) are narrower 
and taller with our model than with ab initio methods, 
which probably reflects the small system size and short 
times of the ab initio simulations compared to ours. In 
summary, our potential reproduces the random tetrahe- 
dral network of amorphous silicon very well, following a 
realistic preparation procedure that starts with a liquid 
phase and cools it down without any artificial changes. 
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VI. CONCLUSION 

We have developed a new empirical potential for sil- 
icon that provides a considerable improvement over ex- 
isting models in describing local structures away from 
equilibrium. The model introduces a theoretically mo- 
tivated functional form that incorporates several coor- 
dination dependent functions to adapt the interactions 
for different coordinations. The fitted potential faith- 
fully reproduces the clastic constants of the equilibrium 
DC structure and also captures the energetics of a wide 
range of point and extended defects and related activa- 
tion energies. The superior description of bulk phases 
and defects is achieved with only thirteen fitting param- 
eters, indicating exceptional transferability of the EDIP 
functional form. Its extended range of success over ex- 
isting models with comparable numbers of parameters 
cannot be attributed to fitting alone. 

For dislocations, this is the first empirical model to 
give a full description of core properties. It predicts the 
correct reconstruction for both the 90°- and 30°-partial 
dislocations, and the reconstruction energies are in agree- 
ment with ab initio data. The bond stretching is also in 
good agreement with ab initio results, pointing to the 
fact that this model predicts reasonably accurately not 
only the energies but also the local structure of disloca- 
tion cores. 

This is the also the first empirical model to predict a 

quench directly from the liquid to the amorphous phase. 
The quality of the resulting amorphous phase, with al- 
most 95% four-fold coordination, is better than with any 
existing empirical preparation method. In some ways the 
amorphous phase with our model is even somewhat closer 
to experiment than with ab initio simulations (surely be- 
cause the latter arc limited to very small system sizes and 
very short times) , which is an encouraging success of the 
empirical approach to materials modeling. 

Our model possesses the same level of efficiency as 
the SW potential, and simulations involving thousands 
of atoms may be readily performed on typical worksta- 
tions. Therefore, it holds promise for successful applica- 
tions to several systems that are still inaccessible to ab 
initio calculations and are outside the range of validity 
of other empirical models. Taking into account the suc- 
cess of our model in describing dislocation properties, we 
suggest that it may also provide a reasonable description 
of small-angle grain boundaries and other such extended 
bulk defects. Considering its success with the amorphous 
and crystalline phases, the model may also describe the 
a-c interface and solid-phase epitaxial growth. 
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TABLE L Values of the parameters that define the potential, 
obtained from a simulated annealing fit to the database described 
in the text. 



A = 7.9821730 eV 
a = 3.1213820 A 
A = 1.4533108 eV 
Qo = 312.1341346 
a = 3.1083847 



B = 1.5075463 A 
c = 2.5609104 A 
7 = 1.1247945 A 
^l = 0.6966326 



p = 1.2085196 
a = 0.5774108 A 
= 0.2523244 
l3 = 0.0070975 



TABLE IL Energy and lattice parameters for high symmetry 
structures. Here we consider the ground-state diamond cubic (DC), 
face centered cubic (FCC), body centered cubic (BCC), simple cu- 
bic (SC) and hexagonal close-packed (HCP) crystals. For DC, the 
cohesive energy per atom E^''^ is given in eV, while for the other 
crystals the difference of the cohesive energy Ec from the ground 
state DC crystal, AE = Ec — Ej^^, is given. All lattice constants ao 
are for the conventional unit cells in A. For the hexagonal crystals 
we also give the c/a ratios. We also compute the lattice constant 
and binding energy of an isolated hexagonal plane (HEX). For this 
comparison we use the SW potential with the isescaled cohesive en- 
ergy for the ground state, as described in Ref.EI. 







DFT/LDA 


EDIP 


SW 


T2 


T3 


DC 


E 


-4.65 


-4.650 


-4.63 


-4.63 


-4.63 






5.43 


5.430 


5.431 


5.431 


5.432 


SC 


AE 


0.348 


0.532 


0.293 


0.343 


0.318 




ao 


2.528 


2.503 


2.612 


2.501 


2.544 


BCC 


AE 


0.525 


1.594 


0.300 


0.644 


0.432 




ao 


3.088 


3.243 


3.245 


3.126 


3.084 


FCC 


AE 


0.566 


1.840 


0.423 


0.548 


0.761 




ao 


3.885 


4.081 


4.147 


3.861 


3.897 


HCP 


AE 


0.552 


0.933 


0.321 


0.551 


0.761 




ao 


2.735 


2.564 


3.647 


2.730 


2.756 




c/a 


1.633 


2.130 


0.884 


1.633 


1.633 


HEX 


AE 


0.774 


0.640 


1.268 








ao 


3.861 


4.018 


4.104 
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TABLE III. Elastic constants Cn, C12, C44, C44 (without inter- 
nal relaxation) and bulk modulus B of the diamond cubic struc- 
ture in Mbar, and the values of two combinations, Cn — C12 and 
C12 — C44, that-are important for stability. The experimental values 
are from Ref.Ej and the tight-binding tesults from Refo. The ab 
initio result for C44 (LDA) is from Ref.LJ. 





EXPT 


LDA 


EDIP 


SW 


T2 


T3 


TB 


Cll 


1.67 




1.75 


1.61 


1.27 


1.43 


1.45 


C12 


0.65 




0.62 


0.82 


0.86 


0.75 


0.85 


C44 


0.81 




0.71 


0.60 


0.10 


0.69 


0.53 


044 




1.11 


1.12 


1.17 


0.92 


1.19 


1.35 


B 


0.99 




0.99 


1.08 


0.98 


0.98 


1.05 


Cll — C12 


1.02 




1.13 


0.79 


0.41 


0.68 


0.60 


C12 — C44 


-0.16 




-0.09 


0.22 


0.76 


0.06 


0.32 



TABLE IV. Ideal (unrelaxed) formation energies of point 

defects (in eV) and relaxation energies AE = E'/"'' - Ey'-"'""^ 
using a 54-a.ijem unit cell. The ab initio (DFT/LP|A) results are 
from Refs.Ej Ej and tight-binding results from Ref.cil. 







DFT/LDA 


EDIP 


SW 


T2 


T3 


TB 


V 


Ef 


3.3-4.3 


3.47 


4.63 


2.83 


4.10 


4.4 




AEf 


0.4-0.6 


0.25 


1.81 


0.02 


0.40 


1.2 


It 


Ef 


3.7-4.8 


6.15 


12.21 


5.85 


6.92 


4.5 




AEf 


0.1-0.2 


2.10 


6.96 


0.82 


3.47 


0.5 


Ih 


Ef 


4.3-5.0 


6.86 


17.10 


5.39 


8.22 


6.3 




AEf 


0.6-1.1 


2.70 


10.15 


1.72 


3.61 


1.3 


CE 


Ef 


5.47 


5.41 


7.90 


6.50 




5.5 




AEf 


0.90 


0.59 


3.26 






1.8 



TABLE V. Unstable Stacking Fault energy (in J/m^) for the un- 
relaxed glide ancl-ejiufile {111} planes as obtained from calculations 
using DFT/LDAE3, SW, and the present interatomic potential. 







DFT/LDA 


EDIP 


SW 


glide 


< 112 > 


2.51 


3.24 


4.78 




< 110 > 


24.71 


13.45 


26.17 


shuffle 


< 110 > 


1.84 


2.16 


1.38 
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TABLE VI. Reconstruction energy (in eV/b) and APD energy (in 
eV) for core structures of partial dislocations, where b is the repeat 
distance of the dislocation. The DFT/LDA,result for reconstruction 
of the 90°-partial dislocation is from Ref.Lj, and for the 30°-partial 
dislocation is from RefEj. Results fot— the interatomic potential 
calculations (SW and T2) are from Re£tj. Tight-binding result for 
the 90°-partial dislocation is from Ref.LJ. 





DFT/LDA 


EDIP 


SW 


T2 


TB 


Reconstruction 












90°-partial 


0.87 


0.84 




0.37 


0.68 


30°-partial 


0.43 


0.33 


0.81 


-0.13 




APD 

90°-partial 
30°-partial 


0.43 


0.41 
0.34 


0.84 


0.37 
-0.13 


1.31 



Direct Quench Artificial Preparation 





LDA 


EDIP 


HM2 


WWW 


ISW 


HMl 


m 


0.2 


0.23 


2.18 


1.2 


0.0 


0.0 


N4 


96.6 


94.43 


93.74 


86.6 


81.0 


73.6 


N5 


3.2 


5.34 


4.04 


11.8 


18.1 


24.6 




0.0 


0.00 


0.01 


0.2 


0.9 


1.5 



TABLE VII. Coordination statistics for model amorphous struc- 
tures at room temperature. A''„ is the percentage of atoms with n 
neighbors closer than the first minimum of g{r). We compare struc- 
tures generated by molecular dynamics simpfetion of a direct quench 
from the liquid using an ab initio methodtj (LDA) and our inter- 
atomic potential (EDIP) with structures generated byrvarious artit 
ical preparation gethods described in the text: IIM2LJ, WWWiIiJ, 
ISWM and HM1L3. 





Pa 


AHa^c 


Z 


ri 




r2 


rs 


e 


0-9 


EDIP 


0.04836 


0.22 


4.054 


2.39 


0.034 


3.84 


5.83 


108.6 


14.0 


EXPT 


0.044-0.054 


< 0.19 


3.90-3.97 


2.34-2.36 


0.07-0.11 


3.84 


5.86 


108.6 


9.4-11.0 


LDA 




0.28 


4.03 


2.38 


0.079 


3.84 




108.3 


15.5 



TABLE VIII. CoBiparison of thermodynamic and strun ti J |ra |l ra gaperties of the present model (EDIP) for a-Si with 
(annealed) ab initi£3 and with (annealed) experimentalEjHLJ'LJ results. Shown are the density pa in A~^, the 
excess enthalpy AHa-c compared to the crystal in eV/atom, the average coordination Z, the mean ri and standard 
deviation of the first neighbor distance in A, the mean second ¥2 and third rs neighbor distances in Aand the 
mean 9 and standard deviation ag of the first neighbor bond angles in degrees. 
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2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 



3.1 3.2 



FIG. 1. The function /(r) that determines the contribution of each neighbor to the 
effective coordination Z. 
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r (A) 

FIG. 2. The two-body interaction V2{r,Z) as a function of separation r shown for co- 
ordinations 3, 4, 6 and 8 and compared with the StilUnger- Weber (SW) pair interaction. 
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FIG. 3. T, which controls the equilibrium angle of the three-body term, as function of 
coordination. 
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FIG. 4. The three-body interaction V;j(r, r, cos Z) for a pair of bonds of fixed length 
r = 2.35 A subtending an angle 0, shown for coordinations 3, 4, 6 and 8 and compared 
with the Stillinger- Weber three-body interaction. 
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FIG. 5. riiKfrgy of the concerted exchange path obtained from calculations using 
DFT/LDAE3'H (diamonds), the SW (dashed line), the Tersoff (dotted line), and the 
present interatomic potential (solid line). The results from ab initio, SW and Tersoff 
are from Ref.LJ. 
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(a) 



(b) 



FIG. 8. Anti-phase defects (APD) for (a) the 30°-partial and (b) the 90°-partial dislo- 
cations. 
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FIG. 9. (a) Pair correlation function and (b) bond angle distribution for the liquid at 
T = 1800 K jfttid zero pressure with the present potential (solid lines) and the ab initio 
model of Ref.E3 (dashed lines). 
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FIG. 10. (a) Pair correlation function and (b) bond angle distribution for the amorphous 
phase at T = 300 K-and zero pressure with the present potential (solid lines) and the ah 
initio model of Ref.EJ (dashed lines). 




r (A) 

FIG. 11. Radial distribution function t(r) — A-nprg{r) for the amorphous phase at room 
temperature and zero pressure using our model, compared with the results of neutron 
scatfetering experiments on pure evaporated-beam-deposited a-Si thin films by Kugler et. 

am. 
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